A Vortex Method for Bi-phasic Fluids 
Interacting with Rigid Bodies 



Mathieu Coquerelle* Jeremie AUard^ Georges-Henri Cottet* 

GRAVIR/IMAG and LMC/IMAG GRAVIR/IMAG LMC/IMAG 

Marie-Paule Cani§ 
GRAVIR/IMAG 



Abstract 

We present an accurate Lagrangian method based on vortex particles, level-sets, and immersed boundary 
methods, for animating the interplay between two fluids and rigid solids. We show that a vortex method 
is a good choice for simulating bi-phase flow, such as liquid and gas, with a good level of realism. Vortex 
particles are localized at the interfaces between the two fluids and within the regions of high turbulence. 
We gain local precision and efficiency from the stable advection permitted by the vorticity formulation. 
Moreover, our numerical method straightforwardly solves the two-way coupling problem between the fluids 
and animated rigid solids. This new approach is validated through numerical comparisons with reference 
experiments from the computational fluid community. We also show that the visually appealing results 
obtained in the CG community can be reproduced with increased efficiency and an easier implementation. 



1 Introduction 

Fluids interacting with solid objects are a common, yet fascinating every-day life experience. Our tendency 
to stare at turbulent liquids and at smoke dynamic behavior, or to observe different objects splashing into 
water makes us very critical when we see such phenomena reproduced with a computer. Yet, the demand 
for plausible fluids in computer-generated movies and games is high. This has led many computer graphics 
researchers to tackle the challenge, although still considered as one of the most difficult problem in the 
computational fluids community. 

The interactions we generally observe in real life involve several components: typically, water with a free 
surface in contact with the air, which interacts with both still and moving rigid bodies. These interactions 
result in really complex phenomena (turbulences, splashes, bubbles, interesting subsequent motion of the 
solids), due to the interplay between the two fluids and the solids. This paper presents a Lagrangian approach, 
based on vortex particles, for accurately, yet efficiently simulating this interplay. We show that vortex particles 
provide a good solution to the modelling of bi-phase flow with a liquid component. Since they concentrate 
the computational power in discontinuous and turbulent regions, the complex phenomena occurring at the 
interface between the two fluids are represented with a good level of precision. Meanwhile, the numerical 
method we present straightforwardly solves the two-way coupling problem between the fluids and animated 
rigid solids, without the need for complex boundary conditions. 
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Previous Work 



The simulation of fluids and of their interplay with solid objects has attracted a lot of attention in the CG 
community within the past few years. Great advances were recently made towards this goal. 

[it) presented the first 3D simulation of liquids. They relied on an Eulerian formulation to solve the Navier- 
Stokes equations for incompressible fluids. A semi-Lagragian method for advecting the fluid which guar- 
anties the stability of the simulation was then proposed by 1 30|. Because of high numerical dissipation, these 
kind of simulations loose part of their vorticity over time, a problem which was tackled in 1 16|. 
Lagrangian methods based on particles provide a good alternative to Eulerian simulations, since they enable 
to dynamically follow the fluid. The Smoothed Particles Hydrodynamics (SPH) formulation was adapted by 
1 12 1 in order to transport an implicit boundary with surface tension. More recently, |26| used this formula- 
tion to simulate fluid-fluid interactions such as boiling water and managed to obtain interesting phenomena 
such as air bubbles in water. Closer to the Navier-Stokes equations for incompressible fluids, 1 28 1 proposed 
a particle based solution called Moving Particle Semi-Implicit (MPS). While straightforwardly able to track 
multiple fluids, particles-only methods mostly suffer from the difficulty to conserve the dynamic properties 
of the fluids, especially when dealing with boundaries. 

Among these approaches, vortex methods 1 18 were used in for animating gaseous phenomena. They 
are particularly interesting since they focus the resolution in the regions of high turbulence. They rely on 
the vorticity formulation of the Navier-Stokes equations, well known in applied mathematics for being an 
accurate alternative to Eulerian methods (see |9|). Vortex particles were also introduced to add subgrid tur- 
bulences on top of an Eulerian simulation of liquids or gazes l29l . impressively counteracting the numerical 
dissipation due to the underlying Eulerian solver. Up to now, no technique was developed in graphics to 
simulate liquids or multi-phases fluids from the vorticity formulation of Navier-Stokes. This may be due to 
the inherent cost of particle methods (from 0{N'^) to 0{NlogN) using (24|), more expensive than Eulerian 
methods since a large number of particles (A^ ^ 1000) are necessary to reach a good level of precision. One of 
the contributions of this paper is to show that when adequately implemented, the vortex formulation achieves 
accurate yet efficient simulations. 

Animation and visualization of water in contact with the air require the precise transport and rendering of 
the interface between them. Because of the discontinuity of the fluids' physical properties at the interface, 
the resolution of the Navier-Stokes equations was most of the time restricted to the liquid component and 
to the interface region [Tl^TTS^. However, this prevents simulating some interesting phenomena such as air 
bubbles inside the liquid. A multi-phase method was recently presented to take care of both fluids, achieving 
impressive bubbles movements |21 1. As in this last work, we simulate bi-phase fluids. We present a different 
solution, based on the vorticity formulation, which solves the problem in an intuitive and efficient way. 

Fluid- structure interaction is considered as a really tough problem by both the CG and applied mathematics 
communities. The most difficult issue is to compute the fluid's velocity at the objects boundaries while 
ensuring that no fluid penetrates any obstacle. 1 19| presented one of the first methods in CG for generating 
interaction forces between fluids and rigid solids. Recently, |5 1 provided a solution that treats the rigid body 
as a fluid and extracts its associated movement from the fluid's velocity. (20 1 presented the first technique 
able to handle the interactions with thin deformable and rigid objects represented by meshes. This method 
was further improved by 1 25 1 to allow the melting and burning of these solids. Closer to our vortex particles 
approach, 1 27 1 used the panel method to impose no-slip and no-through constraints on the fluid by emitting 
new particles at the solid's boundaries. This one-way interaction method, while well adapted for gases, is 
based on an explicit representation of the object which, if the density of particles was too low, does not 
prevent the flow from penetrating. Although achieving impressive results, these methods are sometimes hard 
to implement or require complex treatments where the object touches the fluid. Furthermore, they do not 
provide any validation of the correctness of the fluid's behaviors at the boundaries. Although in the spirit 
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of the rigid fluid method O, the solution we present avoids the expHcit tracking of the fluid/body interface. 
Instead, we use a level set formulation on top of a vorticity creation algorithm to apply forces and account for 
the continuity of velocity at the fluid/body interface. 

Overview 

Our motivation for using vortex particles lies in the following features: they allow to concentrate numerical 
efforts on regions of interest and remain stable for large time-step values. Their main limitation, the associated 
computational cost, can be alleviated by using an auxiliary 3D grid. This is done while keeping the vorticity 
formulation of the equations and remaining physically accurate. 

Our first contribution is to use vortex particles to simulate bi-phase fluids such as liquid and gas. As we show, 
the vortex formulation is very well adapted to tackle this problem, since the vortex particles are created near 
the interface between the two fluids. The second contribution is a novel algorithm to simulate the two-way 
interactions of the bi-phase fluid with animated rigid bodies. Our approach provides a physically sound and 
clear-cut fluid- solid model thanks to an algorithm that is both robust and easy to implement. 

Section previews vortex methods. Section [3l explains how we use them to simulate bi-phase flow, such as 
water in contact with air. Section |4| focuses on the animation of solids interacting with the fluids. Section 
Invalidates the method through numerical comparisons with reference experiments from the computational 
fluid community. Lastly, Section |6l shows that the visually appealing results obtained in the CG community 
can be reproduced with an increased efficiency. 

2 Vortex Particle Methods 

In most fluid simulations, only a part of the flow has an interesting behavior. The vorticity formulation derived 
from the Navier-Stokes equations allows to focus the computational cost on the region of interest using vortex 
particles. 

This section gives a quick overview of the vortex particles method, we refer the reader to 1 9J and the references 
therein for detailed numerical analysis and discussions of this method. 

Definition 

We start from the incompressible viscous 3D Navier-Stokes equations: 

+ (u • V)u - vAu + V/? = (1) 
V-u = (2) 

where is the velocity's time derivative of the fluid's velocity field u, V • u is the divergent of u, p is the 
pressure and v the kinematic viscosity of the fluid. The vorticity is defined as the curl of the velocity: 

CO = V X u 

(note that in the particular case of a 2D fluid, the vorticity is a scalar). Taking the curl of equations and 
^ (after some linear algebra) leads to the vorticity formulation of the Navier-Stokes equations: 

(Ot^{vi-V)(o = ((O- V)u + vAo). (3) 

Solving this single equation is equivalent to solving equations {ij and (|2|). The term {u-V)co represents 
the transport of the vorticity by the fluid's velocity, the left-most term of the right-hand side of equation Q 
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represents the stretching (change of orientation) of the vorticity vector while the latter term represents the 
vorticity diffusion due to viscosity. 

While it would be possible to solve this equation on a traditional 3D grid, the use of particles to transport 
vorticity permits to gain in precision and efficiency. Contrary to velocity which is mostly non zero in the 
whole domain, vorticity is localized in region of turbulences even if there can be high velocities everywhere. 
In consequence, by the use of particles, computational precision is focused where vorticity exists. Another 
advantage is that the advection of vorticity is not subjected to the time- step constraint inherent to Eulerian 
methods, which stands that 5t < \u\tnax- 

Vortex particle methods consist in representing the vector field (0 by a set of particles: 

p 

where Xp, cOp and Vp are respectively the location, strength and volume of particle p and ^ is a smooth 
distribution function, typically a Gaussian. Due to the incompressibility constraint, the volumes Vp remain 
constant. Rewriting equation (|3ll in a Lagrangian formulation, the particles' location and strength are inte- 
grated using: 

-^=u(x„0, (4) 
^ = (a)-V)u + vA(0, (5) 

where Dq/Dt = dq/dt -\- {u- V)q denotes the rate of change of a quantity q in the Lagrangian frame of a fluid 
element advected by the fluid. Particles' velocity and vorticity derivatives have to be determined in a self- 
consistent way from the vorticity field. 



Coupling Vortex Particles with a Grid 

As we plan to simulate fluids which may become turbulent and thus require many particles, we use both 
vortex particles and a 3D grid: firstly, spatial differentiations are cheaper on a grid; secondly it guaranties 
an approximately constant cost when solving the equations; lastly, it solves with no extra cost the problem 
of redistributing the particles over time. This last point is important because particles advected by the fluid 
will naturally tend to cluster in some area or to move away from each other thus leaving unresolved spaces 
between them. 

A fast and accurate solution is to superpose a uniform grid on the particle distribution: at each time step 
particles are remeshed on that grid (thus, the particles' volume Vp = where h is the grid's cells spacing). 
Remeshing at every time-step can be interpreted as a class of high order finite-difference schemes, as long as 
it is done with high order redistribution schemes (see [ Ij for example). 

Given the vorticity field, velocity is computed on the grid by means of a fast grid solver (in our experiments 
we used the FishPack library 1 31 1 in order to solve the following Poisson equation: 

AvA=-(o, (6) 

to obtain the so-called stream function Xj/ which is differentiated to obtain the velocities : 

u = V X VA. (7) 

This solution to compute the velocity field from the vorticity is faster than the particle based solution used in 
1 27 1 when the number of particles increases. 

In practice particles are advanced with a second or fourth order Runge-Kutta method. During each substep 
velocity and vorticity 's time derivative (the right hand side of the vorticity equation Q) are interpolated from 
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the grid onto the particle locations. The accuracy of the overall algorithm heavily relies on the quality of the 
interpolation formulas used to remesh particles on the grid and to interpolate back fields at particle locations. 
It is common practice in CFD to use smooth interpolation formulas which preserve moments of order up to 
2. The stencil on which particles are remeshed extends to the four nearest grid points in each direction. The 
cost of the algorithm thus scales linearly with the number of particles. A number of numerical validations on 
challenging test cases in CFD has allowed to attest the accuracy of remeshed particle methods. These studies 
are backed by a recent theoretical analysis which demonstrate that remeshed particle methods are equivalent 
to a class of high order finite-difference schemes not subject to CFL conditions ((0). 



3 Bi-phase fluids 

We present here an intuitive method for simulating bi-phase fluids which benefits from the vortex particles 
formulation. Before explaining the particularity of bi-phase fluids, we first consider the vorticity problem in 
the presence of a fluid of variable density. 

For variable density and viscosity flows, the vorticity equation Q must be replaced by: 

0)^ + (u • V)(0 = (o) • V)u + V • (vVco) + {Vp X Vp)p~^ 

+ Vx(pg)+V-[(Vx v)Vu] (8) 

where g is the gravity field, p is the density and v = v(p) is the variable kinematic viscosity. This equation 
has to be supplemented by a transport equation for p or for Vp. 

This system is often simplified by assuming small density variations: in the so-called Boussinesq approxima- 
tion, pressure and velocity terms in (|8|) disappear and we are left with 

(0^ + (u-V)(o = ((0-V)u + V-(vV(o)+Vx (pg). (9) 

It is important to note that, in the case of two fluids with constant viscosity and variable density, this equation 
shows that vorticity is produced where density gradients are located, that is at the interface between the 
fluids. Thus computation is localized in this narrow band which clearly reduces the theoretical cost compared 
to traditional methods. Accordingly to the particle's vorticity change in equation ^ becomes: 

= ((0-V)u + V-(vV(o)+Vx (pg). (10) 

We now consider the simpler problem of a bi-phase fluid in a domain £1 where the density (resp. viscosity) 
takes only two different values: pi (resp. Vi) in a domain and p2 (resp. V2) in a domain (^^ = ^^i U ^12)- 
We use a level set (noted 0) to capture the interface between those two domains: 

< for X G ^li , 
(/)(x) { =0forxG^linri2 , 
> for X G ^^2 • 

Typically, we initialize as a signed distance. This level set function satisfies a transport equation which can 
be solved either in its primitive form: 

0, + (u-V)0=O (11) 

or in its gradient (vector) form: 

(V0), + (u • V)V0 = -(V0 • V)u (12) 

The latter equation is very similar to the vorticity equation and thus gradient of (j) and CO are both located near 
the interface. 
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In this case particles carry strengths of vorticity and of from which we can deduce V x p needed in 
equation (ITqI (see below). 

But contrary to 1 8 1, we advect the level set (j) on the grid with a semi-Lagrangian method. The surface tension 
is defined in terms of the curvature of (j) (see 1 6 1): 

TK-C(0)V0. 

where T is the surface tension coefficient. This term is added to the vorticity equation fToll . The level set is 
interpolated onto the particles during the remeshing step. 

In order to be able to solve equation dTol . we need to know the density and viscosity fields. They are defined in 
the whole domain thanks to a smoothed version H of the Heaviside function (with values and 1 respectively 
for positive and negative arugments, and where e is a small parameter of the order of the grid-size) based on 
the level set: 

p = PiH{(^/e) + p2{l-H{(l,/e)). 

The same convention is used for viscosity. From this formulation, it is trivial to deduce the density gradient 
Vp: Vp = (pi - p2)V0(^(0/e)/e. Obtaining V x p is straightforward. 

As the computation of[8lis much more involved and expensive then equation |9j we have implemented the 
latter. Despite the fact that the Boussinesq approximation applies on small density variations, we have found 
that simulating a bi-phase fluid with a high difference of density (e.g. pair ^ 1 and pwater ^ 1000) using this 
technique provides visually realistic results (for both fluid's dynamic and interface localization). 

Algorithm 

In summary, the following process is applied to advance from time tn = nAt to t^^i . Only steps 5 and 8 differ 
from the standard implementation we presented in Section |2| 

1 . Find the stream function by solving equation ^ on the grid, 

2. Compute the velocity u^+^ from equation 0, 

3. Compute the stretching term {co • V)u on the grid, 

4. Interpolate velocity and the stretching term on the particles, 

5. Advance the particles: advect them with velocity (equation ^) and vorticity change (equation fTol)): 
Advect the fluid's level set with the velocity on the grid, 

6. Distribute the particles onto the grid to get the advected vorticity, 

7. Compute and integrate the viscosity term V • (vAco) on the grid, 

8. Create fresh vortex particles, carrying co and V0, from the grid if vorticity is greater than a threshold. 

Note that we used a viscous splitting algorithm 1 7 1 to handle separatly the advection and the diffusion of the 
vorticity. Please refer to figureQfor a visual explanation of the algorithm. 



4 Coupling Fluids with Animated Solids 

We propose a new approach which, with the adjoint power of the vortex particles and level set methods, can 
be used to compute fluid-solid interactions in an intuitive manner and for a low additional cost. Our approach 
is to consider the fluid-body system, during the interaction processing step, as two fluids of different densities 
and different constitutive laws. 

Unlike the method in |5|, our method computing the forces applied by the solid on the fluid (and vice- versa) 
does not require the explicit tracking of the interface. We instead use a level set method combined with a 
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(a) Steps 1 and 2: from vorticity to velocity (on the grid, (b) Steps 4 and 5: interpolation and advection of parti- 
vorticity is denoted by red circles). The magenta curve cles (particles' vorticity is denoted by green discs). The 
represents the interface between the two fluids. thick dashed (resp. thin dotted) magenta represents the 

advected level set (resp. level set at previous time). 
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(c) Step 6 (left figure): particles (in transparent green 
on background) are distributed onto the grid. Step 8: 
remeshing step: create fresh particles where vorticity is 
strong enough (green crosses denote the fact that parti- 
cles were not created because vorticity was too small). 



Figure 1 : Algorithm summary 



vorticity creation penalization term that enforces velocity continuity at the fluid- solid interface. 

In order to simplify the equations we assume here a single solid, the extension to multiple solids being 
straightforward. 

At time zero the body is represented by the zero surface of a level set function noted here (p in order to 
not confuse the reader with the fluid's level set (j). Typically, we initialize the level set function as a signed 
distance to the body boundary, negative inside the solid. This is the only costly part of the method but it does 
not influence its efficiency since it is done as a precomputation. 

We now denote by u"+^ the velocities found in step 2 of the previous algorithm. After this step, we project 
the velocities in the grid onto rigid body velocity U'^+^ inside the solid using the following formula: 

U^+i ^ ( J 5^+1 H{(p/£)dx)/{ J H{(p/£)dx). (13) 

U^^^ stands for the mean velocity inside the solid. A similar equation defines the mean vorticity. Integration 
is performed on the whole domain covered by the solid). At the same time, we enforce continuity of velocities 
at the fluid- solid interface with a penalization technique that we will describe below. 



At this stage, velocity and vorticity in the whole domain are obtained by the formula : 

^^+1 ^ u"+ii^((p/e) +2^+1 (1 -H{(p/£)). (14) 

A similar formula is used to obtain the vorticity field co^^^. Step 3 and 4 of the algorithm now use these fields 
for computing the stretching term and for the interpolations on the particles. 

Step 5 is modified in two ways. First, we take care of the solid's density exactly in the same way as we 
have done for bi-phase fluids. Thus, particles now carry the gradient of a fluid with three different densities. 
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Secondly, we need to advect the solid's level set. This is done just after step 5 of the algorithm of Section 
[3l As the solid's velocity gives a translation and a rotation terms, the latter problem is easily solved by 
simply applying the rigid transformation between the initial level set position at time and its current posi- 
tion. In consequence, this scheme does not suffer from temporal diffusion and only implies a low diffusion 
depending on the order of the spatial interpolation. In our experiments, we use a simple first order interpolant. 

It remains to explain the penalization method to enforce velocity continuity. We use the penalization method 
of |4|: assume we want to solve Navier-Stokes equations in a fluid domain outside a solid domain S, with 
velocity u at the interface. The penalization model reads 

u, + (u-V)u- vAu + V;? = A;^5(a-u), (15) 

where Xs denotes the characteristic function of the solid domain 5" (1 inside and outside) and A 1 is 
a penalization parameter. Adapting this method to the vorticity formulation for our problem leads to the 
substitution of the initial vorticity equation Q by: 

(JL)^ + (u • V) 0) = (o) • V)u + vAco + ?iXs{G) - co) 

+ A5^5nx (a-u), (16) 

where 5^^ is the ID Dirac mass supported by the solid boundary and n is the unit normal, pointed inward. 
Thus, to be consistent with the level set representation for solids, the final equation is: 

(0^ + (u • V)(o = (o) • V)u + vAco + XH{(p/e) {cb - co) 

+ A^^^V(px(a-u). (17) 

The penalization term tends to cancel the vorticity difference inside the body and adds vorticity at the fluid- 
solid interface. This term is computed on the grid after the computation of the fluid's velocity at step 2 of 
the algorithm. It is then interpolated on the particles in step 4. We use an explicit scheme to discretize it and 
choose A = 1 to enforce the stability of this scheme. 



5 Validation 

In this section, we validate the robustness and the physical accuracy of our method for computing the inter- 
actions between fluids and solids. We compare our results to the reference work from |22| obtain with the 
STAR-CD^ software, a well-known software in the Computational Fluid Dynamic community, implementing 
fluid-body interactions from 1 15 1. 

For this purpose, we first study the 2D case of a falling cylinder (which may be seen as a falling disk) in a fluid 
of constant density, submited to gravitational forces. We use the following parameters: boundary conditions 
are periodic in a square domain of size 1.0 x 1.0, cylinder's radius is 0.1, fluid's and cylinder's density are re- 
spectively p fluid = 1 and pcyiinder = 2, kinematic viscosity V is 0.001, gravity is — le^ and £ = 25x = 25y = h 
(£ is used for computing the smoothed Heaviside function and Dirac distribution values). Time steps (resp. 
cells size) are 5tm = 0.01 (resp. hm = 7.8125 • 10"^), 5^256 = 0.0038 (resp. /z256 = 3.90625 • 10"^) and 
5^300 = 0.0027 (resp. hsoo = 3.33 • 10"^) for grids of dimension 128 x 128, 256 x 256 and 300 x 300 respec- 
tively. The two last time steps are constrained by the diffusion stability condition for an exptheslicit scheme. 
Their simulation uses a time step of 5tKem = 0.0005 and a 2D radial grid following the cylinder, 15800 cells, 
80 X 200 nodes (radial x tangential) and first cell at 5r = 0.65 • 10~^. This type of grid is time consuming, 
but well adapted for a sharp resolution of this particular problem as it permits to focus the computations 

^Get more informations on STAR-CD at http://www.cd-adapco.com/sitemap.html 
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Figure 2: y velocities of the cylinder under gravity influence with a RK2 advection for our method, compared 
to the reference curve. 

near the cylinder where precision is the most important. Both ours and Kern's scheme are 2^^ order in space 
and Kern's method is order in time while we are using a Runge-Kutta 2 method for advecting the particles. 

The y component of the velocity obtained with our method (see figure |2|l is quite similar to the one obtained 
with 1 22 1 thus proving that our method well computes the solid's velocity coupled with the application of 
gravity to the solid. As the spatial precision increase, the curves converge to Vy = —0.47. These results 
are quite interesting for computer graphics since, even though we are using time steps 20 times bigger (for 
128 X 128) than the other method and cells 10 times bigger, we still perform a robust treatment of boundaries. 

The vorticity created around the cylinder is also a revealing quantity, both physically (for fluid's dynamics) 
and visually (for turbulences visualisation) important. We compare our vorticity field to the one obtained 
with Kern's method, iso values of (0 are shown in figure |3l vorticity creation and localization are similar. 
Notice that vorticity is not null inside our cylinder, this is because of the small width e used to represent the 
solid. The vorticity trail created behind the cylinder is the signature of a good boundary treatment, we refer 
the reader to 1 23 1 for high-resolution simulations and vorticity analysis of the flow around a cylinder with 
constant velocity. 



6 Implementation and results 

We have implemented our algorithm with a Runge-Kutta 2 advection of vortex particles and second order in 
space differentiations. For the following examples, we have used a low order interpolation for advecting the 
level set. 

In all the results we present now, we use a common gravity g = — lOe^, viscosity for water is v^ater = 
1.0 • 10~^, for air Vair = 0.82 • 10~^, water's density is set to 1 and air's density to 0.001. The surface tension 
coefficient is 10~^. All our results were computed on an Opteron 2GHz with 2Go of memory. 
We present different types of simulations, in figures l^andQcups are falling in water, rendering is done with 
MAYA^^. In figure|6lan initialy vertical volume of water is falling on several solids and in figure|5l we have 
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(a) Reference vorticity isovalues from 1221 . 



(b) Vorticity isovalues using our method with a 
300 X 300 grid. 



Figure 3: Vorticity isovalues at time 0.5 (top) and 2.5 (bottom). 

made two spheres fall in a tank of water, similarly to 1 5 1. For those last two examples, the rendering was done 
with OpenGL. A summary of the performance of the method is provided in table[l] 

In the first example, shown in figure|4| the cup has density pcup = 1.5. We see that, when entering the liquid, 
the cup encapsulates air which makes it roll over after a while, as the captured volume tends to ride up due 
to gravitational forces. One can see in the last two images of the clip, the big bubble escaping and merging 
back with the air. Figure Q shows similar simulations with the same cup starting with different orientations. 
One can observe the effect of surface tension and air capturing on the dynamic of the cups and fluids. The 
grid's dimension is 100 x 100 and the time step used for the two simulations was 0.01. The total time spent 
for 1300 iterations per simulation was approximately 3 hours. 

As a second example, we took a "wall" of water which, under gravity, falls down and breaks the construction 
(see|6ll. One can see in the third image the creation of a wave on the top of the water surface. This wave then 
breaks and merges with the water under it. 

Figure 13 shows two spheres falling in water, this case was inspired from |5 1 and we are using the same grid 
resolution. While obtaining similar dynamics for solids and fluids (the turbulent flow is observable by looking 
at the white "dust" particles), we have computed this sequence with a smaller time step of 0.01, 400 iterations 
were performed in 40 minutes which represent approximately a gain of 58. 
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Figure 4: Interactions of a rigid object with both air and water allows to simulate complex movements. 



Sequence 


Grid Resolution 


Duration 


Time 
Steps 


CPU Time 
/Step 


Cup 1 


100 X 100 X 100 


10 s 


1300 


8.60 s 


Cup 2 


64 X 64 X 64 


10 s 


1000 


2.44 s 


Spheres 


68 X 24 X 292 


4s 


400 


5.96 s 


Pyramid 


80 X 80 X 80 


3s 


600 


12.4 s 



Sequence 


Stream 


Particules 


Grid 


Level- Set 


Rigids 


Rigid 


Surface 


Other 




Solve 








Coupling 


Solver 


Tension 




Cup 1 


32.4 % 


4.93 % 


18.6 % 


23.3 % 


8.28 % 


3.08 % 


5.78 % 


3.66 % 


Cup 2 


32.1 % 


4.60 % 


19.5 % 


21.0% 


8.17 % 


5.38 % 


5.58 % 


3.69 % 


Spheres 


42.2 % 


10.32 % 


21.18 % 


21.2% 


8.35 % 


1.23 % 




3.04 % 


Pyramid 


33,1 % 


2.43 % 


6.81 % 


25.5 % 


24.5 % 


1.06% 


2.05 % 


4.50 % 



Table 1: Performance measurements. Particles stands for particle RK2 advection and distribution, Grid for 
computation on grid (steps 3 and 7), Level-Set for level set advection, diffusion and reinitialiation. 



One can see in the performance synthesis (table [1} that the a third of the computational cost is due to the 
solving of the Poisson equation, while only a small amount of time is dedicated to the particles. The second 
costly part of the algorithm stands for the level set operations. The rigid/fluid coupling takes a little more 
time to compute than grid-based finite differences computations such as the stretching or the viscous term 
computation. 



7 Conclusion and outlook 

We have presented the first vortex particle method that handles a bi-phase fluid interacting with animated rigid 
bodies. Fully based on the vorticity formulation of the Navier- Stokes equations, our method takes benefits of 
the advection of the vortex particles to achieve precision and robustness for large time steps. Meanwhile, it 
relies on a 3D grid for efficiently computing some of the terms and redistributing the particles at each time 
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Figure 5: Two sinking spheres showing drafting, kissing and tumbhng effect 0, reproduced by our method. 

step. Both the bi-phase fluid and the moving soHds are treated using the vorticity equations for a fluid of 
variable density, the associated level sets being advected by the fluid. Moreover, the rigid motion of the solid 
is accurately computed by ensuring the continuity of velocities with the surrounding fluids. 

As our results show, this method brings the two benefits of being a numerically accurate and an efficient way 
to simulate fluids interacting with rigid bodies. 

Our level set approach for modeling the interface between the fluids and immersed bodies could be extended 
to other cases. One can for instance enrich the type of physics underlying the fluid-body interaction, while 
keeping the simplicity provided by the immersed interface approach and its particle discretization. Using 
the framework developed in ITol it is actually possible to handle flexible bodies interacting with fluids, a 
possibility that we plan to explore. Another idea for future work is to implement multi-level particle methods 
in our fluid-body models. One may envision two types of multi-level discretization. One would be to capture 
the interfaces, that is the level- set functions, and the flow, that is the vorticity with different grid resolution. 
Using finer resolution for the level- set function should allow to capture finer details on the interfaces with 
little computational overhead. Indeed it should be emphasized that, since the time-step of particle methods 
is only constrained by the amount of strain in the flow, refining the resolution for the level- set functions does 
not result in a smaller time- step. Alternatively, one may consider using a full (for the flow and the interfaces) 
multi-level vortex method, in the sprit of AMR methods, along the lines of |3 1. Fluid-body interactions is a 
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Figure 6: A pyramid encountering a wall of water, 
challenging problem to clarify the respective benefits of these two approaches. 
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